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Abstract 


In  an  earlier  paper ^  an  algebraic  characterization 
vas  made  of  the  problem  of  resolving  closely  spaced 
plane  waves  incident  on  a  linear  array.  The  character¬ 
ization  suggests  several  data-adaptive  processing  meth¬ 
ods  and  encompasses  the  Wiener,  Maximum  Likelihood,  and 
Pisarenko  methods.  In  this  paper,  the  algebraic  ap¬ 
proach  is  amplified  and  the  results  extended  to  consider 
correlated  noise.  A  recursive  algorithm  is  given  for  a 
particularly  effective  processing  method. 


An  important  array  processing  problem  is  that  of 
determining  the  directions  of  propagation  of  plane  waves 
incident  on  a  linear  array  of  uniformly  spaced  sensors. 2 
Contemporary  spectral  analysis  has  led  to  the  develop¬ 
ment  of  several  array  processing  methods  that  are  able 
to  resolve  plane  waves  with  nearly  identical  directions 
of  propagation.  These  methods  include  the  Wiener  Pre¬ 
diction  Filter  method, 3  the  Maximum  Likelihood  method, 3 
and  the  Pisarenko  method.1*  This  paper  amplifies  and 
extends  an  algebraic  approach  given  earlier1  based  upon 
an  algebraic  characterization  of  the  array  processing 
problem.  The  results  encompass  the  methods  mentioned 
above  and  include  the  case  of  correlated  noise.  A  re¬ 
cursive  algorithm  is  presented  for  implementation  of  a 
particularly  effective  processing  method. 

Model  of  the  Array  Data 

Consider  the  complex  sinsoidal  time-space  plane 
wave  T{t,r)  as  represented  by 

f(t,r)  =  Ae  -  U) 

where  A  is  the  complex  amplitude,  t  is  the  continuous 
time  variable,  r  *  xi+yj+zk  is  the  continuous  space^ 
variable,  (v  is  the  (temporal)  frequency,  and  k^  »  kxi 
+k  J+k  k  is  the  wavenumber  (spatial  frequency).  This 
wave  travels  in  the  direction  of  -k  with  a  speed  of 

propagation  c  =  yjy  .  Let  us  now  monitor  this  wave  with 

a  linear  array  placed  along  the  x-axis  whereby  y  ■  z  ■  0 
as  shown  in  Figure  1.  The  detected  signal  is 


J[wt+k  x] 

fx(t,x)  *  Ae  .  (2) 

From  this  ideal  data  it  is  theoretically  possible  to 
determine  the  values  of  the  parameters  u  and  kx-  Fur¬ 
thermore,  if  the  speed  of  propagation  is  a  known  con¬ 
stant  or  a  known  function  of  frequency,  we  can  then 
calculate  the  wavenumber's  magnitude  from 

W  •  (3) 

Because  we  do  not  have  com[ lete  knowledge  of  the  wave¬ 
number  it ,  we  cannot  unambiguously  determine  the  direc¬ 
tion  of  propagation.  However,  we  can  determine  the 
polar  angle  y  associated  with  the  wave  as 


This  angle  defines  a  cone  whose  central  axis  lies  along 
the  linear  array.  This  information  alone  is  sufficient 
for  many  applications.  For  example,  the  wave  may  be 
known  a  priori  to  be  traveling  in  the  xy  plane.  This  is 
the  case  we  shall  consider  hereaftecsr  - _ ;  .  — 
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Figure  1 .  Linear  array  alotl 


We  have  reduced  the  problem  from  four  dimensions  ir, 
the  variables  to  two  dimensions,  through  the  constraint 
of  a  linear  array.  We  can  further  reduce  the  problem  to 
one  dimension  by  noting  that  time  and  space  are  inde¬ 
pendent  quantities  in  this  model.  Thus  we  can  perform 
our  analyses  for  u  and  kx  separately.  Ordinary-  time 
series  processing  such  as  filtering  or  spectral  estima¬ 
tion  can  be  applied  to  each  sensor  to  first  determine 
the  presence  of  signals  at  a  particular  temporal  fre¬ 
quency  u.  These  outputs,  one  for  each  sensor,  are  then 
spatially  processed  to  determine  the  directions  of 
sources  radiating  at  the  frequence  u.  Hereafter,  we 
shall  suppress  the  time  domain  and  consider  only  the 
3  5 

spatial  dimension. 

Figure  2  depicts  a  linear  array  of  p  sensors  uni¬ 
formly  spaced  d  units  apart.  A  plane  wave  is  impinging 
upon  the  array  with  an  incident  angle  8.  Hoting  that 
the  incident  angle  is  complementary  to  the  polar  angle, 
we  have 

k  k  k 

sin  8  -  T-2r  =  «  TT  1 

|k|  w/c 

from  which  it  is  seen  that 

k  .  ^-B-fn-6  (5) 

x  A 

where  A  is  the  wavelength  of  the  plane  wave.  Defining 
our  origin  at  sensor  zero,  the  nbh  sensor  will  sample 


origin  a 
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the  wave  at  the  point  x  »  nd.  Hence,  at  any  particular 
instant  in  tine  the  array  output  ia,  from  (2) 


yin) 

where  t  ia  a 

in:'*  .'lilt . 


■  Ae  L  *  J  (6) 

o*n«p 

phase  angle  dependent  on  the  sampling 
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Figure  2.  Plane  wave  with  incident  angle  0. 


"tie  set  of  p  instantaneous  spatial  samples  (6)  as 
measured  by  the  array  is  referred  to  as  a  "snapshot." 

Ir.  this  case  the  snapshot  is  a  sampled  complex  sinusoid 
whore  sampled  spatial  frequency  is  given  by 

-  ■  ?:J  ?in  9  ,  Clearly,  an  estimate  of  the  sinusoid's 
frequency  directly  yields  an  estimate  of  the  direction 
of  rrepagation  of  the  plane  wave  if  the  wavelength  X  is 
Kr.  wu.  As  such,  spectral  estimation  is  seen  to  play  a 
imminent  role  in  linear  array  processing. 

We  now  generalize  our  model  to  include  multiple 
p  i  -it. •  waves  incident  on  an  array  in  which  the  sensor 
in  i  i  ."it  ions  are  contaminated  by  white  measurement  noise. 
I  f  t  here  are  a  total  of  q  plane  waves  and  the  k**1  plane 
wave  ha.  a  direction  of  propagation  0^,  it  follows  by 
.'up  rpnr.ition  that  the  snapshot  will  have  the  form 


\  "a  Jhu. 

y(n)  «  n(n)  +  y  A^e  e  ,  0£n<_p-l 

k=l 

where  the  q  cinuosoid  frequencies  are  given  by 
2nd  sin  0. 


(7) 


and  n(n)  are  uncorrelated  zero  mean  random  variables 

n 

with  variance  o‘  that  represent  the  measurement  noise. 

We  assume  that  the  w,  are  all  different. 

Our  objective  is  to  estimate  the  frequency  param¬ 
eters  ui  using  these  snapshot  measurements.  We  are 
particularly  interested  in  the  ability  to  resolve,  or 
distinguish  between  two  plane  waves  with  very  similar 
frequencies  (i.e.,  w  i  w^).  This  estimation  capability 
require.-,  the  utilization  of  a  number  of  snapshots  taken 
sequentially  in  time.  Our  data  then  has  the  form 


(  n )  =  n  ( 


n_ln) 
Pi 


A  e 
k 


J*km  Jn<1,k 


1  ^m  £M,  0  <_n  <_  p-1 


k=l 


(8) 


where  m  ir  tin-  snapshot  index  and  M  i3  the  total  number 
of  snapshots.  In  this  model,  we  assume  that  the  phase 


angles  are  uncorrelated  random  variables  uniformly 
distributed  on  [-*,  +  *].  Their  random  nature  arises 
from  the  Independence  of  the  sinusoidal  sources  and  from 
the  approximate  randomness  of  tlae-sampllng  far  below 
the  Nyqulst  rate. 

It  la  convenient  at  this  point  to  represent  the 
problem  in  vector  notation.  We  represent  the  mth  snap¬ 
shot  (8)  by  the  p*l  column  vector 

y  *  ty  (0)  y  (1)  ...  y  (p-1))'  (9) 

"in  m  in  rn 

and  define  the  pure  complex  sinusoid  vector  as 

S  -  [1  ej“  eJ?“  ...  eJ(p'1)“]'  .  (10) 

*Tl) 

Lastly,  the  noise  vector  associated  with  the  mth  snap¬ 
shot  is  defined  as 

th  “  tB(l)  •••  hB(p-l)]'  •  (11) 

With  the  above  notation,  we  may  compactly  represent  the 
snapshots  by  the  data  vector  equation 


*■  "  An  *  2 
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k-1 


(12) 


The  array  data  (12)  is  random  due  to  its  dependency 
on  the  random  phase  angles  and  the  contaminating 
noise  n^n).  Assuming  that  these  random  variables  are 
pairwise  uncorrelated  and  statistically  invariant  with 
respect  to  the  snapshot  index  m,  it  follows  that  each 
data  vector  v  is  a  vindoved  realization  of  a  vide-sense 
stationary  raKdom  process.  The  mean  value  of  this  proc¬ 
ess  is  the  zero  vector. 


UzJ  *  E^}  + 
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vhile  its  covariance  matrix  is  specified  by 


(13) 
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where  1^  is  the  p*p  identity  matrix  and  *  lA^!' 

the  power  of  the  kth  incident,  plane  wave.  Since  th, 
random  vector  process  is  vide-sense  stationary,  the 
covariance  matrix  R  must  be  positive  ser.i-definite, 
Toeplitz,  and  Hermitian. 

We  now  describe  three  contemporary  array  process  ir..- 
methods  and  then  present  an  algebraic  approach  to  iden¬ 
tifying  the  frequencies  (w^),  based  upon  the  structure 
of  the  data  ^  and  the  associated  covariance  natrix  R. 

Contemporary  Processing  Methods 


Wiener  Filter  Method 


The  Wiener  Filter  method  is  bar.rd  on  filtering  the 
data  such  that  the  signal-to-noise  ratio  at  the  filter 
output  is  maximized. 8  it  is  essentially  a  linear  pre¬ 
diction  approach  that  is  quite  sir.ilar^to  the  Maxima- 
Entropy  method  of  spectral  estimation.'  Many  adaptive 
array  processing  algorithms  are  equivalent  to  the  Wiener 
Filter  method,  including  Alan's  orthonom.il  lattice 
filter  algorithm.^  For  the  array  processing  problem  an 
optimum  weighting  vector  is  obtained  by*' 


a  =  uK'-'h  (lb) 

wliere  h*[l  0  0  ...  0 1  ’  and  u  is  nil  arbitrary 

scalar.  As  in  the  Maximum  gntropv  Method,  the  power 
spectrum  may  be  computed  by 


pw(“) 


S+a  a+S 


(16) 


Maximum  Likelihood  Method 


The  Maximum  Likelihood  method  is  based  on  filtering 
the  data  such  that  power  at  the  frequency  of  interest! a 
passed  and  all  other  frequency  components  are  rejected 
in  an  optimal  manner.  In  our  notation,  the  power  spec¬ 
trum  is  given  by  6*9,10 


PMLW 


s+s_1s 
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(17) 


Pisarenko  Method 


The  Pisarenko  method  has  not  found  as  widespread 
use  as  the  Wiener  Filter  method  and  the  Maximum  Likeli¬ 
hood  method.  Haykin"  recently  applied  part  of  the 
Pisarenko  method- 1  to  the  array  processing  problem  via 
a  special  autoregressive-moving  average  (ARMA)  model. 

The  complete^!  sarenko  method  is  based  on  a  theorem  of 
Caratheodory  that  allows  decomposition  of  the  exact 
truncated  covariance  sequence  r(n),  O^n^q-l,  into  a 
positively-weighted  sum  of  q  complex  sinusoids  and  white 
noise.  The  method  has  three  steps: 

(i)  identifying  and  removing  the  noise  contri¬ 
bution  to  the  covariance  matrix, 

(ii)  forming  the  qxq  covariance  matrix  Rq  and 

and  analysing  the  single  eigenvector  corre- 
pondlng  to  the  unique  minimum  eigenvalue 

1  .  *  0  to  determine  the  sinusoid  frequencies, 

min 

and 

(iii)  solving  a  set  of  q  simultaneous  linear 
equations  for  the  sinusoid  powers. 

Algebraic  Processing  Approach 

We  now  formulate  a  generalized  minimization  problem 
which  suggests  several  particular  methods.  Under  dif¬ 
ferent  constraints,  the  solution  of  this  problem  encom¬ 
passes  each  of  the  methods  in  the  previous  section. 

Let  us  consider  a  general  nontrivial  p*l  coeffi¬ 
cient  vector  a  that  is  orthogonal  to  the  noise-free 
component  of  each  of  the  data  vectors  From  (12), 

this  orthogonality  is  defined  by  the  inner  product 
relationship 

0  * 

■J]V  <■•  8,  >  U8) 

k«l  k 

Since  the  (u.  }  are  all  different  and  the  (♦,  )  are  ran- 
K  Kin 

dom  in  nature,  a  little  thought  will  convince  oneself 

that  a  must  be  orthogonal  to  each  of  the  q  sinusoid 
vectors  S  ,  1 <  k <  q.  We  next  define  the  general  z- 
^k 

transform  A(z)  of  the  coefficient  vector  a  by 

A(z)  »  <a,  z»>  (19) 

where  z  *  [l  z*1  z_?  ...  z1_p]  .  It  is  then  readily 

shown  that  the  orthogonality  of  a  to  each  S  ,  1  <  k  <  q, 

”**k 

implies  that  A(z)  must  have  q  finite  zeros  located  on 

J“k 

the  unit  circle  at  the  points  z^  =  e  ,  l^k^q.  With 
this  in  mind,  the  required  frequencies  can  be  deter¬ 
mined  by  examination  of  the  zeros  of  A(z). 


In  the  Idealized  noise-free  case,  the  snapshot  vec¬ 
tors  are  members  of  the  q-dimensional  subspace  span¬ 
ned  by  the  q  linearly  independent  sinusoid  vectors  S  , 

”®k 

l£k£q.  Thus,  forq<pthere  always  exists  a  p*l  vector 
a  that  is  orthogonal  to  the  noise- free  components  of  the 
snapshot  vectors  y^.  If  q>_p,  there  generally  does  not 
exist  such  a  vector  a.  Furthermore,  when  noise  is 
present,  even  if  q<p,  there  generally  doe:  not  exist  a 
vector  a  that  is  orthogonal  to  each  of  the  noise- 
contaminated  snapshot  vectors  jfa.  nonetheless,  it  is 
intuitively  desirable  to  select  a  coefficient  vector 
which  is  nearly  orthogonal  to  each  of  the  snaps hot  vec¬ 
tors  in  some  well-defined  manner,  and  to  determine  the 
plane  wave  sinusoidal  frequencies  by  examination  of  the 
zeros  of  the  z-transform  of  this  coefficient  vector.  A 
convenient  method  to  evaluate  these  zero  locations  is  to 
search  for  nulls  in  the  magnitude  of  the  Fourier  trans¬ 
form  of  the  coefficient  vector.  Since  there  can  be  more 
zeros  than  plane  waves  (i.e.,  p-l>q),  we  can  estimate 
the  (P^)  in  order  to  separate  "signal  zeros"  from  "noise 
zeros". 

To  obtain  a  mathematical  measure  of  closeness  to 
orthogonality,  it  is  beneficial  to  define  an  orthogo¬ 
nality  error  vector  e(a)  whose  mth  element  is  the  inner 
product  of  a  and  y^.  denoted  by  <  a ,yr >  .  We  define  the 
optimum  a  to  be  the  vector  a°  that  minimizes  some  posi¬ 
tive  definite  functional  f  of  £(a).  Hence  we  write 

e(a)  =  [e(l)  e(2)  ...  e(M)]1 

where 

e(m)  =  <  a,^  > 

and 

f(e(a°))  *  min  f(e(a)) 
atA 

where  A  is  a  constraint  set  for  the  solution  vector  a°. 

A  constraint  is  generally  necessary  for  the  minimization 
to  be  well-defined,  that  is,  for  a°  to  be  unique  and 
nontrivial. 

We  must  choose  an  inner  product  for  (20)  and  ar. 
error  functional  f  for  (21).  Let  us  choose  in  particu¬ 
lar  the  standard  vector  inner  product  <  a,^  >  *  a'y  in 
which  case  we  have 


(20) 

(21) 


(22) 


A  convenient  positive  definite  functional  for  an  error 
vector  is  the  mean  square  error  criterion 

f(e)  •  E{  ||  e  || P)  (23) 

where  ||  e  ||  »  ^|e(l)|P  ♦  ...  ♦  |e(M)|P,  the  Euclidean 

norm  of  e.l**  It  will  be  computationally  expedient  to 
normalize  this  criterion  by  the  length  M  of  the  vector 
e.  From  (22)  and  (23)  we  have 


■  a*Rn  (2l) 
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where  R  is  the  covariance  matrix  defined  in  (lit).  Apply¬ 
ing  (21),  the  quadratic  form  (2it)  must  now  be  minimized 
according  to  some  constraint  that  causes  t»°  to  be  unique 
and  nontrivial.  Next  we  consider  two  possible  con¬ 
straints,  a  linear  constraint  and  a  quadratic  con¬ 
straint. 

Linear  Constraint 

The  first  constraint  is  that  a0  lies  on  a  hyper- 
planc  specified  by 

A  •  (at  CP:a+h+hta  »  2)  (25) 

where  h  is  a  nontrivial  pxl  vector  that  characterizes 
the  orientation  of  the  hyperplane.  The  solution  to 
(21)  with  this  constraint  is  given  by 


1°  *r-rVlR'1ii 

h  R  h 


and  the  criterion's  minimum  value  is 

f(a°)  «  *  .  (27) 

h  R  h 

Quadratic  Constraint 

The  second  constraint  is  that  a°  lies  on  the  quad¬ 
ratic  surface  specified  by 

A  *  (a <  CP:a+Wa  -  1)  (28) 

where  W  is  a  positive  definite,  Hermitian  matrix  which 
characterizes  the  quadratic  surface.  The  solution  to 
(21)  with  this  constraint  is  given  by 


/ x+ .  Wx  . 
i  -min  -min 


and  the  criterion's  minimum  value  is 


Extension  to  Correlated  Noise 

Through  selection  of  the  constraint  set,  the  alge¬ 
braic  approach  extends  quite  readily  to  the  case  in 
which  the  contaminating  noise  is  correlated.  Such  is 
the  case  when  the  noise  is  due  not  only  to  sensor  meas¬ 
urement  noise  but  also  to  a  directional  background 
noise  field  in  the  array  environment.  Note  that  any 
undesired  signal  (jasning  interference,  for  example) 
may  be  considered  as  correlated  noise. 

Generalizing  (lL),  we  have  that  the  data  covariance 
matrix  for  the  case  of  correlated  noise  is  given  by 


R  *  o  B  ♦ 


s+ 

-v-k 


where  the  noise  covariance  matrix  B  is  defined  by 

o2B  •  E(n  r/)  .  (32) 

TBmtK 

We  assume  that  the  shape  of  the  noise  spectrum  is  known, 
which  implies  knowledge  of  B. 

Returning  to  the  quadratic  constraint  (29,30),  we 

note  ( X  .  ,  x  .  )  is  the  solution  to 
min*  -min 

(R-X  .  W)x  ,  *  0  .  (33) 

min  -min  — 

We  have  seen  that  for  the  choice  W  =  I  we  have  the 

P 

Pisarenko  solution.  In  this  case  it  is  well  known  that 

2 

X  ,  =o  and  that  we  essentially  have  a  white-noise 
min  -« 

power  cancellation  algorithm.  Hence  it  is  a  simple 
step  to  choose 


.  .  and  achieve  a  colored-noise  correlation  cancellation 

algorithm.  This  step  can  be  Justified  further  by- 


algorithm.  This  s 
rewriting  (31)  as 


(R-o2b) 


S  S+ 
^V“k 


where  (X  ,x  ,  )  is  the  minimum-eigenvalue  and  eigen- 
min  -min 

vector  pair  of  W^F. 

The  above  solutions  encompass  the  three  contempo¬ 
rary  methods  described  earlier.  For h  *  [l  0  ...  0]', 
Z;uation  (26)  is  precisely  the  Wiener  Filter  solution 


h+R-1h 


Just  as  in  linear  prediction,  this 


constraint  implies  that  the  first  element  of  a0  is 
fixed  at  1  and  the  other  elements  are  unconstrained. 

For  h  =  F  ,  (27)  is  precisely  the  Maximum  Likelihood 
solution.  This  constraint  requires  A°(z)  to  have  unity 
gain  at  z  *  eJ“,  while  the  minimization  strategy  opti¬ 
mally  reduces  the  gain  at  other  frequencies.  For 
W  =  Ip,  the  quadratic  surface  is  the  hypersphere  of 
radius  one,  and  (29)  is  a  generalized  version  of  the 
Pisarenko  method.  There  are  several  differences. 

First,  no  special  ARMA  model  is  invoked,  as  is  done  by 
Haykin.1'  Second,  neither  noise  power  removal  nor  ma¬ 
trix  order  reduction  are  required,  as  they  are  In  the 
Pisarenko  method.  Third,  this  method  is  based  upon  a 
minimization  strategy  and  so  Justifies  estimates,  gen¬ 
erally  even  non-Toeplitz ,  of  the  covariance  matrix  R. 

In  the  special  case  of  a  Toeplitz  estimate  matrix,  a 
power  identification  technique  similar  to  the  Pisarenko 
method  ear  be  employed,  as  is  shown  later.  Fourth,  the 
general  constraint  matrix  W  allows  greater  flexibility 
in  the  solution.  The  quadratic  constraint  solution 
generalizes  the  Pisarenko  method  and  extends  it  to  the 
multi pie- snapshot  array  processing  problem. 


and  remembering  that  we  seek  a  vector  that  is  orthogonal 

to  the  sinusoid  vectors  (S  ).  Also,  it  is  apparent 
'“k 

+  * 

from  the  constraint  a  Ba  =  1  that  we  are  again  specify¬ 
ing  a  set  of  vectors  with  constant  norm,  but  now  the 
norm  is  determined  by  the  noise  covariance  matrix  B. 

The  linear  constraint  (26,2 7)  may  also  be  extended 
for  correlated  noise.  We  shall  only  consider  the 
Wiener  solution,  since  it  has  been  shown  to  achieve  _ 
greater  resolution  than  the  Maximum  Likel  ihood  Method.  ’ 
We  will  show  later  that  to  extend  the  Wiener  linear  pre¬ 
diction  solution,  a  reasonable  constraint  h  given  B  is 


Note  that  we  no  longer  have  a  linear  predictor. 

Example 

To  illustrate  the  linear  and  quadratic  constraint 
solutions,  let  us  consider  the  case  of  a  single  plane 
wave  of  power  F.  and  spatial  frequency  u  incident  on 
an  array  of  two  sensors  in  a  correlated  noise  field. 
For  this  case  the  noise  covariance  matrix  is  given  by 


2  2 
o  B  =  a 


[::] 


where  |b|<l,  and  the  data  vector  covariance  matrix  by 


2 

Pe  1  «•  o  b* 


'^e  1  ♦  oT) 


First  we  consider  the  linear  and  quadratic  con¬ 
straint  solutions  without  accounting  for  the  correlated 
noise  (i.e.,  the  Wiener  and  Pisarenito  methods).  For 
the  linear  constraint  we  choose  =  [l  0  ...  0]' 

and  find  from  (26)  that 


P  e  1  +  o^b 


Pefir.irg  the  signal -to-noise  ratio  by  SHR 
see  that  A°(z)  has  a  zero  located  at 


L I  SNR  +  b  e  1 

L  SKR  +  1  J 


>1  2  , 

IP^  1  ♦  oZb| 


Thus  A°(z)  has  a  zero  located  at 


OT  e  +  b 

Z  *  J 

Ju. 


We  aee  in  this  case  that  the  zero  lies  directly  on  the 
unit  circle,  regardless  of  the  signal-to-noise  ratio. 
In  fact,  when  the  noise  is  white,  the  zero  perfectly 
indicates  the  plane  wave  spatial  frequency  u^.  How¬ 
ever,  when  the  noise  is  correlated,  there  is  a  fre¬ 
quency  bias  present  that  again  becomes  greater  as  the 
signal-to-noise  ratio  decreases. 

Let  us  consider  now  the  linear  and  quadratic  con¬ 
straint  solutions  which  account  for  the  correlated 
noise.  For  the  linear  constraint  we  choose 
h  »  B(l  0  ...  0)'  and  find  that 


Ve  A-b) 

P^l-be  1)  +  o(l-|b| 


where  ii  is  a  scalar  function  of  P.  ,  u  ,  o  ,  and  b. 
Thus  A°(z)  has  a  zero  located  at  1  1 


.1  SNPtl-be  *) 

-Jw.  p 

SHR(l-b  e  A)  +  l-|b| 


As  in  (1.0),  we  see  that  pure  white  noise  will  still 
cause  the  zero  to  migrate  away  from  the  unit  circle, 
and  that  correlated  noise  will  introduce  frequency  bias. 
However,  as  the  noise  becomes  "more  correlated"  (i.e., 
|b|  ■* 1 ),  the  zero  moves  closer  to  the  unit  circle  and 
asymptotically  indicates  the  exact  plane  wave  spatial 
frequency  regardless  of  the  signal-to-noise  ratio. 
Hote  that  the  effect  of  an  interfering  harmonic  source 
J“>p 

(i.e.,  b  •  e  )  is  completely  removed. 

For  the  quadratic  constraint  we  choos  V  =  B  arid 
find  that 


As  has  been  noted  elsewhere, *5  this  linear- 
prediction  solution  suffers  from  zero  migration  away 
from  the  unit  circle  even  when  the  noise  is  white  (i.e., 
b  =  0).  This  migration  degrades  resolution,  since 
applying  the  Fourier  transform  to  evaluate  zero  loca¬ 
tions  may  indicate  only  a  single  null  when  in  fact 
there  are  two  zeros  close  together  somewhat  off  the 
unit  circle.  Furthermore,  we  see  there  is  a  frequency 
bias  introduced  by  the  correlated  noise.  This  bias 
becomes  greater  as  the  signal-to-noise  ratio  decreases. 

For  the  quadratic  constraint  we  choose  W  =  I  and 
find  'Yom  (29)  that  ** 


J“l  2 
P,e  +  o  b 


where  v  Is  a  Bcalar  function  of  and  b.  Thus  A°(z) 
has  a  zero  located  at 


We  see  that  the  zero  indicates  the  exact  plane  wave 
spatial  frequency  «  ,  regardless  of  the  signal-to-noise 
ratio  or  the  particular  value  of  b.  For  this  reason, 
we  expect  the  quadratic-constraint  solution  to  obtain 
high  resolution. 

To  summarize  the  development  to  this  point,  the 
algebraic  approach  is  based  on  approximating  an  orthog¬ 
onality  condition  between  a  solution  vector  and  each 
of  the  data  vectors.  This  approach  encompasses  three 
contemporary  array  processing  methods  and  readily 
extends  to  the  case  of  correlated  noise. 


Implementation  of  the 


dratic-Constralnt  Solution 


In  the  previous  section  we  saw  that  the  quadratic- 
constraint  solution  (29,30)  is  a  promising  array  proc¬ 
essing  method  in  terms  of  its  perfect  resolution  given 
exact  covariance  values.  However,  it  requires  an 
eigenvalue-eigenvector  computation  that  seems  to  be 
quite  burdensome.  Fortunately,  a  simple  recursive  algo¬ 
rithm  can  be  derived  using,  the  nature  of  the  array 
processing  problem.  , 

First  we  recall  the  standard  "inverse  iteration"11' 
method  for  finding  the  minimum  eigenvalue  and  eigen¬ 
vector  pair  of  a  complex  matrix  D.  Consider  the  se¬ 
quence  of  vectors  (x^)  defined  by 

D*w  ‘  \-i  (1,T) 

where  Xq  is  a  nonzero  and  arbitrary.  Ac  K  increases, 
we  have 


5 


It  can  be  shown  that  6X  Is  approximately  given  by 


^k  Tain 


and 


Vk-1 

t 


*min 


This  method  is  appropriate  to  the  array  processing  prob¬ 
lem,  in  which  the  data  arrives  sequentially.  Assume 
that  from  M  snapshots  we  have  estimated  the  covariance 

matrix  by  R^  and  obtained  the  desired  pair  (^,^1, 

When  the  next  snapshot  is  available,  we  form  RM+^  and 
compute  from  (1*7)  using  D  *  RM+1  and  x^  »  x^. 

Since  the  inverse  iteration  method  generally  has  fast 
convergence,  a  single  Iteration  of  (1*7)  for  may  t>e 

sufficient  as  long  as  is  only  slowly  time-varying. 

To  accelerate  convergence,  we  can  apply  the  "inverse 
iteration  of  Wielandt"1^  wherein  an  approximation  of 
*min  is  subtracted  from  the  main  diagonal  of  D  before 

iterating.  Given  R„  . ,  we  use  X  to  approximate  Xu  . 

M+l  m  1 

The  iteration  is  given  by 


(RM+1  ■  XMI)5M+1  *  ’ 

3m*1  ■*  *nin  *  and 


^I+A+l 


=  A  -*■  A  .  . 

M+l  min 


(1*8) 


For  a  Toeplitz  and  Hermitian  covariance  matrix  estimate, 

2 

each  iteration  can  be  performed  with  0( p  )  multiply- 
adis  using  Zohar's  algorithm.17  An  alternate  algorithm 
by  Oueguen1®  can  be  used  to  avoid  numerical  difficulties 


that 


may  be  associated  with  (RM+1  -  *MI). 

For  the  case  of  correlated  noise,  we  have 


D 


r-i: 


'  M+l  ’ 


iteration  to 
iteration  is 


We  suggest  generalizing  the  accelerated 

avoid  calculation  of  P  Our  general 
given  by 


(RM+1  ~  V  'b+l  ^  =  % 


M  ’ 


4l+l  "  4un  > 

M  t 


M*1 


and 


min 


(1*9) 


.  x  (4D)x 

«>  ■  - - -  (r?) 

xWx 

Our  approximation  to  the  new  X  is  then  given  by  the 

sum  of  4X  and  the  previous  X  ,  .  With  appropriate 

Bln 

definitions,  this  approximation  replaces  X  in  (1*9) 
when  we  expect  the  new  covariance  matrix  estimate  to 
differ  considerably  from  the  previous  estimate. 

Relationship  to  Linear-Constraint  Solution 

The  iterative  implementation  of  the  quadratic- 
constraint  solution  gives  insight  into  the  linear- 
constraint  solution.  Namely,  the  first  iteration  of 
(1*7)  with  D»R  and  x^  *  h  yields  the  linear-constraint 
solution  of  (26)  within  a  constant  of  proportionality. 
Repeated  calculation  of  the  linear-constraint  solution, 
with  h  at  each  step  equal  to  a°  of  the  previous  step, 
is  in  fact  an  iterative  implementation  of  the  quadratic- 
constraint  solution.  It  is  apparent  that  at  each  step, 
the  constraining  hyperplane  is  realigned  according  to 
the  estimated  solution.  With  these  Insights,  it  is 
reasonable  to  choose 


(53) 


as  the  linear  constraint  for  correlated  noise,  since  it 
yields  the  first  step  of  the  iterative  quadratic- 
constraint  solution  (without  acceleration)  for  corre¬ 
lated  noise  given  in  (1*9).  This  Justifies  the  choice 
made  earlier  in  (36). 

In  this  section  we  have  presented  a  recursive 
algorithm  (1*9)  for  implementing  the  quadratic-constraint 
solution.  The  algorithm  makes  use  of  the  sequential 
nature  of  the  snapshot  data  to  efficiently  employ  in¬ 
verse  Iteration.  The  algorithm  includes  the  case  of 
correlated  noise.  A  modification  to  the  algorithm  (52) 
was  presented  for  the  case  where  successive  covariance 
matrix  estimates  differ  considerably. 

Covariance  Matrix  Estimate 


2 

Each  iteration  still  only  requires  0(p  )  multiply-adds, 

and  B  is  never  calculated. 

In  some  applications,  it  may  not  be  desirable  to 
calculate  the  solution  vector  after  every  snapshot. 

For  instance,  forming  a  new  covariance  matrix  estimate 
and  calculating  a  new  solution  only  every  L  snapshots 
reduces  the  average  computation  rate  by  a  factor  of  L. 
Unfortunately,  if  the  new  covariance  matrix  estimate 
differs  considerably  from  the  previous,  the  previous 
eigenvalue  may  be  a  poor  approximation  to  *mjn*  and 

convergence  will  be  slowed.  To  obtain  a  better  eigen¬ 
value  approximation,  we  apply  perturbation  techniques.^ 
Suppose  that  0  and  W  are  Hermitian  matrices  and  that  we 
have  solved  the  eigenvalue-eigenvector  problem 

Dx  *  XWx  .  (50) 

Applying  a  Hermitian  perturbation  40  to  D  we  have  the 
new  problem 

(D+4D)(x+4x)  *  (X  +  4X  )W(x  ♦  4x)  .  (51) 


To  employ  the  proposed  processing  methods,  an 
estimate  of  the  covariance  matrix  is  required.  From 
this  estimate,  a  solution  vector  i  ;  obtained  and  the 
zeros  of  the  vector's  z-transform  examined  to  determine 
the  plane  wave  spatial  frequencies.  Given  a  p*l  solu¬ 
tion  vector  and  q  plane  waves,  q<p,  there  will  be  q 
"signal"  zeros  and  p-q-1  "noise"  zeros.  These  zeros 
must  be  separated  from  one  another.  It  is  well  known 
that  in  the  linear  prediction  solution,  dominant  fre¬ 
quency  components  will  generate  zeros  closer  to  the  unit 
circle  than  less  powerful  components;  thus,  a  simple  way 
to  evaluate  signal  zero  locations  is  to  search  for  nulls 
in  the  solution  vector's  Fourier  transform.  For  the 
quadratic  solution  in  white  noise,  it  can  be  shown^^ 
that  all  of  the  zeros  will  be  on  the  unit  circle  when 
the  covariance  matrix  estimate  is  both  Hermitian  and 
Toeplitz.  Thus  the  estimated  frequencies  can  be  di¬ 
rectly  employed  in  a  power  determination  technique'1  *  ' 
and  the  zeros  separated  on  a  basis  of  signal  power  as 
before . 
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3 

A  standard  covariance  matrix  estimate  is 
M 

vsEyl-  (53) 

m=l 

This  estimate  is  unbiased.  Hermit ion,  but  in  general 
not  Toeplitz.  Furthermore,  only  one  lag  product  from 
each  data  vector  is  used  in  formulating  each  element  of 
RM-  An  alternate  estimate  is  the  matrix  R^  whose  ele¬ 
ments  are  given  by 

RM(i,j)  «=  c(i-j),  l<i,j<p  (5M 

where 

M  p-n-1 

c(n)  *  m]C^  2  ya(t+n)yI(tK  °i"-p-1 

m=l  1=0 

c(n)  =  c#(-n)  ,  -p  +  l±n<0. 

Tills  estimate  is  unbiased,  Hermitian,  and  Toeplitz. 
Also,  p-n  lag  products  from  each  data  vector  are  used 
in  formulating  each  element  c(n).  Thus  the  estimate  of 
(51*)  has  lower  variance  than  that  of  (53). 

In  the  following  simulations,  the  standard  non- 
Toeplitz  estimate  1^(53)  will  be  used  in  the  linear- 
constraint  solution  in  order  to  compare  with  previous 
simulations. 3  However,  for  the  quadratic-constraint 
solution  the  Toeplitz  structure  is  important,  hence  the 

Toeplitz  estimate  R  ( 51* )  will  be  used. 

M 

Simulation  Results 

To  compare  the  performance  of  these  two  processing 
methods,  the  data  vectors  (12)  were  generated  by  com¬ 
puter  simulation.  The  simulation  model  corresponded  to 
that  chosen  by  Gabriel^  in  his  comparative  paper. 

Namely,  the  case  of  two  sources  incident  on  an  array 
with  white  noise  was  considered.  The  parameter  selec¬ 
tions  were  q=  2,  p  =  8,  a2  =  1 ,  =  A?  =  31-62  (30  dB  SNR) 

and  3.162  (10  dB  SNR),  6^18°,  62  =  22°,  M=50  (many 
snapshots)  and  10  (few  snapshots),  and  d  =  X/2 . 

With  vhite  noise,  the  linear  solution  and  the 
quadratic  solution  correspond  to  the  Wiener  and 
Pisarenko  methods.  The  simulation  results  are  shown  in 
Figure  3.  In  this  figure,  the  linear  solution  has  been 
evaluated  via  its  Fourier  transform  and  the  quadratic 
solution  via  the  power  determination  technique.  Over- 
layed  solutions  for  ten  different  realizations  of  the 
random  data  are  shown  to  give  a  sense  of  each  method's 
consistency. 

These  results  show  that  both  methods  work  well  at 
high  SNR  with  many  snapshots.  However,  the  linear 
solution  performs  poorly  at  low  SNR  with  few  snapshots, 
while  the  quadratic  solution  continues  to  give  good 
resolution  and  good  suppression  of  spurious  effects. 

In  general,  the  quadratic  solution  has  shovn  better 
performance  than  the  linear  solution  over  a  vide  range 
of  conditions. ^  Further  simulations  are  underway  to 
compare  the  performance  of  the  methods  in  correlated 
noise  and  to  evaluate  the  recursive  algorithm  presented 
above.  The  results  will  be  given  at  the  conference. 

Conclusions 

We  have  detailed  an  algebraic  approach  to  array 
processing  based  upon  approximation  of  an  orthogonality 
condition.  The  approach  encompasses  several  contempo¬ 
rary,  high  resolution  methods.  Previous  results  were 
extended  to  the  case  of  correlated  noise,  and  a  recur¬ 
sive  algorithm  presented  for  the  quadratic-constraint 


solution.  The  quadratic-constraint  solution  appears  to 
he  particularly  effective  and  suggests  further  investi¬ 
gation  of  eigen-analysis  array  processing  methods  and 
their  implementation. 
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Figure  3.  Two-source  simulation  with  sources  at  18  and  22  degrees. 

(A)  30  dB  SNR,  50  snapshots  | 

(B)  10  dB  SNR,  10  snapshots  /  Linear  solution 

(C)  Single  trial  from  (B)  ) 

(0)  30  dB  SNR,  50  snapshots  } 

(E)  10  dB  SNR,  10  snapshots 

(F)  Single  trial  from  (E) 
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developed.  This  characterization  describes  many  well-known  techniques 
for  solving  this  problem  as  well  as  Introducing  a  new  effective  approach. 
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